clear all
use "Stones.dta"

set type double

***********************************
* Drop countries with obs less than 160 out of 180
***********************************
tabulate c, generate(ccc)
forvalues i = 1/107 {
di "ccc"`i'
count if kg4 !=. & ccc`i' ==1
drop if r(N) <= 160 & ccc`i' ==1
}
drop ccc*
***********************************

* Rename
gen kjpy = kyen4 // Stones = 4 
gen kg = kg4
* Share
egen csum_kjpy = sum(kjpy), by(t) // country-sum
gen sh = kjpy/csum_kjpy
gen lnsh = ln(kjpy/csum_kjpy)
* Value
gen lnva = ln(kjpy)
* Price (unit value)
gen lnpr = ln(kjpy/kg)
* Exchange Rate
gen lner = -ln(erlccjpy)
* Panel Data
xtset c t
* Dummy variables
tabulate t, generate(t)
tabulate c, generate(c)

***********************************
* Benchmark Regression (IV)
***********************************
xtivreg2 lnsh t1-t180 (lnpr = l.lner), fe endog(lnpr)
nlcom (sigma: 1-_b[lnpr] )
***********************************
* Benchmark Regression (LS)
***********************************
xtreg lnsh t1-t180 lnpr, fe
nlcom (sigma: 1-_b[lnpr] )


**********************************
* Seach the point of normalization 
* (Optional)
**********************************
/*
forvalues tt = 1/180{
quietly gen byte baseyear = 1 if t == `tt'
by c (baseyear), sort: gen lnsh1 = lnsh[1]
gsort t c
quietly reg lner lnsh1 if baseyear == 1
*scalar tval = _b[lnsh1]/_se[lnsh1]
scalar rss = e(rss)
scalar obs = e(N)
**********************************
scalar omegahat = -_b[lnsh1]
scalar tauhat = - _b[_cons]
predict del, residuals
by c (baseyear), sort: gen del1 = del[1]
gsort t c
drop del 
gen lnuh = tauhat + omegahat*lnsh // - del1
gen lnph = lner + lnuh // <========
quietly xtivreg2 lnsh t1-t180 (lnph = lner), fe endog(lnph)
*di `tt' " " tval " " obs " " e(cdf)
di `tt' " " rss " " obs " " e(cdf)
*nlcom (sigma: 1-_b[lnph])
drop baseyear lnsh1 del1 lnuh lnph
}
*/
**********************************

bysort c(t): center lnsh lner lnpr  // horizonstal differences from time-average
bysort t(c): center c_lnsh c_lner c_lnpr // vertical differeces from item-average

**********************************
*SUR
**********************************
quietly gen byte baseyear = 1 if t == 7 // t = Theta
by c (baseyear), sort: gen lnsh1 = lnsh[1]
gsort t c

reg lner lnsh1 if baseyear == 1 // REGSUPPLY  OMEGA
estimate store sup

reg c_c_lnsh c_c_lner // KAPPA
estimate store demex

suest sup demex
nlcom (1 - [demex_mean]_b[c_c_lner]/(1-[sup_mean]_b[lnsh1]*[demex_mean]_b[c_c_lner]) )

**********************************
* FEENSTRA
**********************************
gen yy = (c_c_lner)^2
gen xx = (c_c_lnsh)^2
gen xy = (c_c_lner)*(c_c_lnsh)

ivregress gmm yy (xx xy = c1-c19), nocons
nlcom (sigma1: (_b[xy]+sqrt(_b[xy]^2+4*_b[xx]))/(2*_b[xx])+1    )
*nlcom (sigma2: (_b[xy]-sqrt(_b[xy]^2+4*_b[xx]))/(2*_b[xx])+1    )

scalar rho1 = (_b[xy]+sqrt(_b[xy]^2+4*_b[xx]))/(2)
scalar rho2 = (_b[xy]-sqrt(_b[xy]^2+4*_b[xx]))/(2)

gen nu1 =c_c_lner - c_c_lnsh*rho1
gen lnnu1 = ln(nu1)
gen expnu1 = exp(nu1)
ivreg2 c_c_lnsh (c_c_lner = expnu1 lnnu1), endog(c_c_lner)
*ivreg2 c_c_lnsh (c_c_lner = nu1), endog(c_c_lner)
nlcom (sigma: 1-_b[c_c_lner])
